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'nI" ■ Abstract 

We present numerical verification of hyperbolic nature for chaotic attractor in 

Q ■ a system of two coupled non-autonomous van der Pol oscillators (Kuznetsov, Phys. 

U: Rev. Lett., 95, 14410f , 2005). At certain parameter values, in the four-dimensional 

^ I phase space of the Poincare map a toroidal domain (a direct product of a circle and 

r— H . a three-dimensional ball) is determined, which is mapped into itself and contains 

1 the attractor we analyze. In accordance with the computations, in this absorbing 
I domain the conditions of hyperbolicity are valid, which are formulated in terms of 

^ ■ contracting and expanding cones in the tangent spaces (the vector spaces of the 

Tj- , small state perturbations). 

O ; 

^ . Mathematical theory of chaotic dynamics based on a rigorous axiomatic foundation 

O ■ exploits a concept of hyperbolicity [1-8]. 

^ I An orbit in phase space of a dynamical system is called hyperbolic if there are tra- 

' jectories approaching exponentially the original orbit, and those departing from it in a 

• ^ ■ similar manner. Moreover, an arbitrary small perturbation of a state on the original or- 
' bit must admit representation via a linear combination of the growing and the decaying 
> ! perturbations. 

I In dissipative systems contracting the space volume the attractors may occur, which 

^ I consist exclusively of the hyperbolic orbits. These are attractors with strong chaotic prop- 
erties, like existence of the well-defined invariant SRB-measure, a possibility of description 
in terms of Markov partitions and symbolic dynamics, positive metric and topological 
entropy etc. Such hyperbolic (or, more definitely, uniformly hyperbolic) attractors are 
robust or structurally stable, that means insensitivity of the type of dynamics and of the 
phase space structure in respect to slight variations of functions and parameters in the 
evolutionary equations. 

Although the basic statements of the hyperbolic theory were formulated 40 years 
ago, no convincing examples of physical systems were introduced with uniform hyper- 
bolic attractors. In textbooks and reviews on nonlinear dynamics, such attractors are 
represented by artificial mathematical constructions, like Plykin attractor and Smale - 
Williams solenoid [1-8]. For realistic systems, in which the chaotic dynamics is math- 
ematically proved, like the Lorenz model [9,10], the strange attractors do not relate to 
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the class of uniformly hyperbolic (not all axiomatic statements of the classic hyperbolic 
theory are valid for them). Some aspects of possible existence of hyperbolic attractors in 
differential equations were discussed e.g. in Refs. [11-14]. 

In a recent paper of one of the authors [15], an idea was advanced of implementation of 
a hyperbolic attractor in a system of two coupled non-autonomous van der Pol oscillators. 
In a Poincare map that determines evolution on a period of the external driving, a chaotic 
attractor has been found, which demonstrates some characteristic signs of hyperbolic 
attractors. By a nature of transformation of the phase space volume in a course of the 
evolution over a period, it is similar to the Smale - Williams solenoid. It looks robust: 
the Cantor-like transverse structure and the positive Lyapunov exponent are insensitive 
to variation of parameters in the equations. An analogous system has been built as an 
electronic device and studied in experiment [16]. 

Obviously, it would be desirable to have a mathematical confirmation of the hyper- 
bohc nature of the attractor. As Sinai has suggested in due time [1] , one possible way for 
substantiation of the hyperbolicity for attractor of a Poincare map consists in numerical 
verification of certain sufficient conditions formulated in terms of inclusion for expanding 
and contracting cones in tangent vector space (the space of small perturbation vectors). In 
this paper, we discuss a method and present results of computer verification of these con- 
ditions in application to the chaotic attractor in a system of two coupled non-autonomous 
van der Pol oscillators. 

The system proposed in Ref. [15] is represented by a set of differential equations 



It consists of two subsystems, the van der Pol oscillators with characteristic frequencies 
ujq and 2ujq. Here x and u represent coordinate and velocity for the first oscillator, and y 
and V for the second one. In each oscillator the parameter responsible for the birth of the 
limit cycle, is forced to swing slowly with period T and amplitude A. As the parameter 
modulation is of opposite phase, the subsystems generate turn by turn, each on its own 
half-period T. The coupling is characterized by parameter e. The first oscillator affects 
the second one via a quadratic term in the equation. The backward coupling is introduced 
by a product of the variable y and an auxiliary signal of frequency u^o- It is assumed that 
the interval T contains an integer number of periods of the auxiliary signal A^^o = ujoT/2tt, 
so the external driving is periodic. For a detailed study, we select 



Qualitatively, the system (P) operates as follows. Let the first oscillator on a stage 
of generation have some phase ip: x cc sin(co'o^ + V")- The squared value contains the 
second harmonic: cos(2u;o^ + S'?/'), and its phase is 2iIj. As the half-period comes to the 
end, the term effects as priming for the excitation of the second oscillator, and the 
oscillations of y get the phase 2?/'. Half a period later, the mixture of these oscillations 
with the auxiliary signal stimulates excitation of the first oscillator, which accepts this 
phase 2ip. Obviously, on subsequent periods the phase of the first oscillator will follow 
approximately the relation 



X = uqu, u = —uqx + {Acos2nt/T — x'^)u + (e/cuo)?/ coscuot, 
ij = 2uqv, V = —2uQy + {—Acos2nt/T — y^)v + {e/2uJo)x^. 



(1) 



ujo = 2tt, T = 6, A = 5, e = 0.5. 



(2) 



ipn+i = 2V^n + const (mod27r). 



(3) 
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(Here the constant accounts a phase shift in a course of transfer of the excitation from one 
oscillator to another and back.) The relation (jS)) called the Bernoulli map is well known 
as one of the simplest model examples in the chaos theory.^ 

For accurate description of the discrete time dynamics, we turn to the Poincare map 
[2-8, 17,18]. Let us have a vector x„ = f as a state of the system 

at tn=nT. From solution of the differential equations with the initial condition x„, we 
get a new vector x„+i at tn+i = {n+l)T. As it is determined uniquely by x„, we introduce 
a function that maps the four- dimensional space {x,u,y,v} into itself: x„+i = T(x„). 

This Poincare map appears due to evolution determined by differential equations with 
smooth and bounded right-hand parts in a finite domain of variables {x, u, y, v}. In accor- 
dance with theorems of existence, uniqueness, continuity, and differentiability of solutions 
of differential equations, the map T is a diffeomorphism, a one-to-one differentiable map 
of class [17]. 

Further, we will deal always with description of the dynamics in terms of the Poincare 
map. In particular, under the phase space we mean the four-dimensional space {x, u, y, v}, 
with X, u, y, V relating to an instant t„. An orbit means a discrete sequence of points in 
this space; attractor is an invariant attractive set composed of such orbits etc. 

In a course of iterations of the map x„+i = T(x„), we have expansion of a small phase- 
space volume in a direction associated with the phase in the approximate equation Q 
and contraction in the rest three directions. Interpreting the mapping geometrically, let 
us imagine a solid toroid embedded in the 4-dimensional space (a direct product of a circle 
and a three-dimensional ball) and associate one iteration of the map with longitudinal 
stretch of the toroid, with contraction in the transversal directions, and insertion of the 
doubly folded "tube" into the original toroid. It is analogous to the construction of Smale 
and Williams with the only difference that we deal with four-dimensional rather than the 
three-dimensional phase space. 

The mentioned toroid will be referred to as an absorbing domain U. It means that 
under application of the map T the images of all points from U belong to its interior: 
T(f/) C Int U. The attractor may be defined as intersection of the images of the original 

oo 

domain under multiple action of the map: A = T"(f/). 

n=l 

To write down an analytic expression for the domain U it is convenient to redefine the 
coordinate system. We introduce new variables {xq, Xi, X2, x^} as follows: 

^0 = 3;/ro, Xi = (u — Cuxx)/ri, X2 = y — Cy^x — CyuU, x^ = v — Cy^x ~ CvuU — Cyyy. (4) 

To determine the constants we accumulate a large number of points {x, u, y, v} on the 
attractor in the Poincare section by numerical solution of the equations (^. Then, by 
the least square method we find out the coefficients to minimize the mean-square values 
< {u - Cuxx)"^ >, < iy - CyxX - CyuuY > , < [v - CyxX " c^„M - CyyyY >. Gcometrically, 
it corresponds to directing the coordinate axes along the principal axes of ellipsoid that 
approximates the attractor. Additionally, we normalize Xq and Xi by appropriate factors 

^The constant in Eq. Q may be removed by a shift of origin for the phase variable. We stress that the 
phase ip cannot be defined globally on the whole time interval T: it has sense only in the context of the 
discrete time description. Indeed, on the stage when the first oscillator does not generate, its amplitude 
is small, and phase is not well defined. 
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to have ]xli=]x'li^l/2. Finally, at the parameter set we get 



■vu 



■ux 



0.438, c, 
0.029, c, 



■vy 



yx 



-0.042, = 0.226, c^^ = -0.218, 
-0.118, ro = 0.812, ri = 0.721. 



(5) 



In the new coordinates, let us define the absorbing domain U by the inequality: 



Empirically selected constants in this expression are r = 0.94, dr = 0.4, d = 0.15. Figure 
1 gives evidence that this is indeed an absorbing domain. For initial points distributed 
over a border of U we perform numerical solution of the differential equations on an 
interval T and plot the results in the coordinates 



As the whole figure is placed inside the unit circle Rf + R2 = 1, the images of the initial 
points belong to the interior of U. 

In Fig. 2 we show a three-dimensional projection to illustrate mutual location of the 
domains U and T{U). It is analogous to that considered on the first step of the construc- 
tion of the Smale - Williams attractor: take a torus ("a plastic doughnut"), stretch it 
twice, contract transversally, fold twice and squeeze into its original volume. The trans- 
formed "doughnut" T{U) looks like a narrow band because of very strong compression of 
the phase volume in respective directions in a course of the evolution. 

We will verify hyperbolicity conditions required by a theorem (see e.g. [1, 13]) adopted 
for the problem under consideration. Unlike the general formulation, it is sufficient for 
us to deal with a diffeomorphism of class C°° in the Euclidian space M.^ {xq, Xi, X2, x^}. 
That is the Poincare map T(x). Evolution of a perturbed state x -|- (5x corresponds to 
transformation of the perturbation vector 5x in linear approximation 5x' = DTx5x, where 
DTx is the Jacobi matrix at x: DTa, = {dx'Jdxj}, i,j = 0, 1, 2, 3. The notion DT~^ 
designates the derivative matrix for the inverse mapping T~^(x). 

Theorem [1,13]. Suppose that a diffeomorphism T of class C°° maps a bounded 
domain U C R'^ into itself: T{U) C Int U, and A clntU is an invariant subset for the 
diffeomorphism. The set A will be uniformly hyperbolic if there exists a constant 7 > 1 
and the following conditions hold: 

1. For each x G A in the space V^; of 4D vectors 5x the expanding and contracting 
cones and may be defined, such that ||DTxU|| > 7 ||u|| for all u G S"^, and 
||DTx"'^v|| > 7 ||v|| for all v G C^; moreover, for all x G A they satisfy S*^ fl = 
and SZ + C^^ = V^. 

2. The cones 5*^ are invariant in respect to action of DT, and are invariant in respect 



to action of BT ^ i.e. for all x G A DTx(52) C S^,. and DT^^(Q) c C^_i 



If the formulated conditions are valid for a whole absorbing domain containing the 
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attractor, say. 



T"(f/), they are obviously true for the attractor A = T"(f/). 



n=l 
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Let us consider in some detail the procedure of computer verification of these condi- 
tions. Been given a point x = {a;o, Xi, X2, x^} G U, we perform numerical solution of 
Eqs. Q on the interval t G [0, T] with the initial state 

x\t=o = f^oXo, u\t=o = riXi + CuxX, 

y\t=Q = X2 -\- CyxX + CyuU, v\t=o = X3 + CyxX + C^uU + C^yU 

and get the image 

X = T(x) = {xq, Xi, 3^25-^3} 

/tq) i^u CuxX )/ri, y CyxX CyuU , v c^xX c^uU c^yy }. 

(9) 

In parallel, we solve numerically the linearized equations for vectors of small perturbations 
over the same period. In the original variables they are 

6x = uJo6u, 6u = —LJo6x — 2xu6x + {Acos27it/T — x'^)Su + {e/ujo)6ycosujot, , , 
5y = 2uJov, 6v = —2uJo6y — 2yv5y + {—Acos2-7Tt/T — y'^)5v + {£/uJo)x6x. 

Passage to the redefined coordinates and back may be done with the relations 

6x0 = 6x/ro, 6x1 = {6u - CuxSx)/ri, 6x2 = 6y - CyxSx - Cyu6u, 
6x3 = 6v — CvxSx — c^u^u — Cvy6y. 

The equations pO|) are solved along the orbit started at x for four times, each time 
with such an initial vector u= {6xi} that unity is posed in a row from to 3, and other 
elements are zero. Then, we get four vector-columns and compose a matrix U = DTx of 
them. 

Starting at x, an initial perturbation vector u after one iteration of the Poincare 
map yields u' = Uu. A squared Euclidean norm of this vector is ||u'||^ = u-^U^Uu, 
where T means the transposition. Using the inverse matrix U^^ we can write as well 
u = U~^u' and ||u||^ = u'"^U~^'^U~^u'. A condition that u' represents an image of 
a vector belonging to the expanding cone 5*^, is given by an inequality ||u'|| > 7||u||, 
or u'^ (U^^'^U^-^ ^ 7~^) ^' < 0- Starting at x' = T(x) = {xq, x[, X2,x'^}, an initial 
vector u' transforms into u" = U'u', and we have ||u"||^ = u'^U'^U'u'. The ex- 
panding cone >S':^(x) ^'^ = T(x) is determined by an inequality ||u"|| > 7||u'||, or 

u'^ (jJ'^V - 72) u' > 0. 

The 4x4 matrix U'"^U' is positive definite and symmetric. Let dg, di, d2, be vectors 
of the orthonormal basis. The matrix D = (do, di, d2, da) transforms the matrix U'^U' 
to diagonal form: 

D^U'^^U'D == {A^^,,} (12) 

Let the eigenvalues be enumerated in decreasing order. As we have one expanding and 
three contracting directions, then, Aq > 1 and Af 23 < 1- Now, we suppose that 7 is 
selected in such way that Aq > 7^ and Af 2.3 < 7^-^ Then, in the matrix 

D^(U'^U'-7')D = {(A^-7')5,,} (13) 



^This property is checked in a course of computations at each analyzed point of the absorbing domain 
naturaUy: its violation would entail a non-correct operation of taking a square root of a negative number. 
The inequalities for eigenvalues of the matrix UjUx ensure fulfillment of the condition that a sum of 
subsets of the linear vector space (that is a set of all possible linear combinations of vectors from the 
expanding and contracting cones) is the full 4D vector space: + = V. 
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we have one positive and three negative elements on the diagonal. By an additional scale 
change along the coordinate axes 

S = {s;%j}, So = V^Ag-y, .1,2,3 = v^^^^ (14) 
it is reduced to the canonical form 

S^D^(U'^U' - f )DS = H' = {h!,5ij}, K = h K^,^, = -1. (15) 
A vector c={l, 01,02,03} belongs to the expanding cone S^^^y if c^H'c > 0, or 

cl + cl + cl< 1. (16) 

In the 3D space {01,02,03} it corresponds to interior of the unit ball. 

With the same transformations, the matrix U~^'^U~^ — 7"^ takes a form 

S^D^(U-i'^U-^ - 7-2)DS = H = {hij}. (17) 

(Note that it is symmetric: hij = hji.) A vector €={1,01,02,03} represents an image of a 
vector belonging to the expanding cone, if c^Hc < 0, or 

^00 + ^01 Cl + /l02C2 + /iosCs + hiQCi + hiic\ + /li2Ci02 + /ll3Ci03+ , . 

/l20C2 + h2lCiC2 + /l22ci + /l23C2C3 + /IsqCs + /^aiCiOs + /l32C2C3 + /issOg < 0. 

In the space {01,02,03} it corresponds to interior of some eUipsoid. The inclusion DTx(5'2) C 
'^T(x) n^sans that the ellipsoid has to be placed inside the unit ball. To formulate a suffi- 
cient condition for this, we determine a center of the ellipsoid from the equations 

hiiCi + /li2C2 + /ll3C3 = -/llO, 

/l2lCl + /l22C2 + /l23C3 = -/l20, (19) 
^3lCl + /l32C2 + /l33C3 = -/i30, 

and estimate a distance of this point from the center of the ball: 



4 + 4 + 4- (20) 

With a transfer of the origin to the center of the ellipsoid, the equation for its surface 
becomes 

hiic\ + hi2CiC2 + /il3ClC3 + /l2lClC2 + /i22C2 + /i23C2C3 + hiClCs + /l32C2C3 + /I33C3 = , (21) 

where Cj = Oj — q, and 



= -(/loo + hoiCi + ho2C2 + ho3C3 + hioCi + hii4 + hi2CiC2 + /ll3Cl03+ 
/l20C2 + h2lCiC2 + h224 + /«'23C2C3 + HsqCs + hsiCiCs + /l32C2C3 + hsacl). 



(22) 



Now, we consider a symmetric 3x3 matrix h — {hij}, i, j = 1,2,3. In the diagonal 
representation of this matrix, under appropriate orthogonal coordinate transformation 
(ci, 02, C3) — > (^1, ^2, C3)) the equation of the ellipsoid surface becomes 

+ 12(^ + 13^^1^, (23) 
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where /i, I2, h are eigenvalues of h. The largest semiaxis of this ellipsoid is expressed via 
the minimal eigenvalue: 

A sufficient condition for the ellipsoid to be positioned inside the ball is given by an 
inequality 

r-max + P < 1. (25) 

It completes the procedure of verification of the expanding cones inclusion for the point 

X. 

It may be shown that with 7 <1 the application of the above procedure in U is 
equivalent to verification of the condition in the domain xGT^(f/) for contracting cones 
with the parameter 7' = I/7 > 1: 'DT~^{Cx'^) C C^-i^^^y It is so because the cones S"^ 
and C^/'^ are complimentary sets: S"' U C^^'^ = V. (Here S'^ with 7 < 1 corresponds to the 
cone of vectors, which either expand, or contract, but no stronger than by the factor 7.) 
Hence, fulfillment of the inequality ()25|) checked inside U for two parameters 7 and I/7 
would imply that both conditions for expanding and for contracting cones are valid in the 
domain T^(f/), which contains the attractor.^ This is sufficient to draw a conclusion on 
the hyperbolic nature of the attractor. 

The computer verification of the required inclusions for the expanding and contract- 
ing cones was performed at the parameter values Q in the coordinate system Q, 
Computations of the Poincare map and of the Jacobi matrices were produced by means 
of joint numerical solution of the differential equations together with linearized equa- 
tions (fnH) on the time interval T. We used the Runge - Kutta method of the 8-th order 
based on formulas of Dormand and Prince with automatic selection of step (the accuracy 
for one step was assigned to be 10~^^) and an extrapolation method (the accuracy for 
one step assigned 10^^^) [19]. For solution of sets of linear algebraic equations, matrix 
diagonalization, and eigenvalue problem solving, we used sub-programs from the library 
LAPACK [20]. 

In accordance with our computations, at 7^ = 1.1 the sufficient condition (j^^ of 
correct inclusion for the expanding cones DTx(S'2) C S'^i^^-^ is valid in the whole absorbing 
domain U. To discuss details, let us consider a 3D hypersurface defined by an equation 

2 

+ {x2/df + (xs/df = R". (26) 

At i? = 1 it corresponds to a border of the domain [/; at i? < 1 it belongs to its interior. 
We can parametrize this hypersurface by three angle coordinates 0, ip, and 6: 

xq = {Rdr cos 9 + r) sin ip, xi = {Rdr cos 9 + r) cos ip, /r,7\ 
X2 = Rd sin 9 cos (f), X3 = Rd sin 9 sin 0. 

The variable may be regarded as a phase of the first oscillator at the Poincare cross- 
section, and as a phase of the second oscillator at the same instant. Numerical computa- 
tions on a 3D grid with step 271 /M at M = 50 show that the value rmax + P = f{R, 0, V'; 
at fixed R depends essentially on ip and 9, while dependence on is very weak. On the 

^At the same 7, the cones S"' and have a common border only at 7 = 1, while at 7 > 1 they do 
not intersect, as required by the theorem conditions: S^CiC^ = 0. 







/dr 
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plot of the function / one global maximum can be seen of value varied in dependence 
on (j) and R. At i? = 1 and some (f) the maximum reaches fmax ~0. 929441 (that corre- 
sponds to a point M on the border of the domain U with coordinates Xq = —0.102628, 
Xi = —0.544957, X2 = 0.000581, 0:3 = 0.040066), but remains definitely less than 1, see 
Fig. 3^ Panel (b) illustrates mutual disposition for the cones DTx(5'2) and 5'^(x) the 
point M. The plot shows a 3D cross-section of the 4D vector space Vt(x) by a hyperplane 
orthogonal to the expanding direction. The coordinate axes are principal semiaxes of the 
ellipsoid representing the cross-section of the cone •S':^^^)- scale selection along the 

axes, it looks like a ball. The ellipsoid representing the cross-section of DTx(S'^) looks like 
a narrow "needle" , because of high degree of phase volume compression in two directions. 
Its disposition inside the large ball testifies the condition DTx(S'2) C S'^^^y The ball 
circumscribed around the ellipsoid is posed inside the large ball too; that expresses the 
sufficient condition (j25|) . For smaller R the global maximum of rmax + P only decreases 
(Fig. 4a). Analogous computations with other values of 7 indicate that the required in- 
clusions for the cones S take place at least in the interval 0.64 < 7^ < 1.35 (Fig. 4b). As 
explained, the correctness of the condition with 7 < 1 implies the condition for the con- 
tracting cones DT-^(Cx^^) C C}^-^^) in the domain T'^{U). We conclude that in T^{U), 
both conditions for expanding and contracting cones are true, say, at 7^=1.1.^ Hence, 
the analyzed attractor is uniformly hyperbolic. This assertion, although not proven in 
a classic mathematical style, follows with definiteness from the theorem, conditions of 
which have been checked in the computations. Assuming the hyperbolicity established, 
let us illustrate now some attributes of the hyperbolic dynamics. 

To start, we note that dynamics on the attractor is chaotic. In a course of time 
evolution, both oscillators generate turn-by-turn, passing the excitation one to another. 
Figure 5 shows typical plots for x and y obtained from numerical solution of Eqs. (P). 
Panel (a) presents a single sample, and panel (b) shows five superimposed samples of 
the same signal on successive time intervals. Fig. (b) gives evidence that the process is 
not periodic. In fact, it is chaos, which manifests itself in irregular displacement of the 
maxima and minima of the waveforms relative to the envelope on successive time intervals 
T. 

To have a quantitative indicator of chaos we turn to Lyapunov exponents. With mul- 
tiple iterations of the Poincare map and Jacobi matrix computations, we trace evolution 
of four perturbation vectors by means of their subsequent multiplication by the Jacobi 
matrices obtained in a course of the evolution. At each iteration, the Gram - Schmidt 
orthogonalization and normalization are performed for the set of vectors. The Lyapunov 
exponents are determined as the mean rates for growth or decrease of the accumulating 
sums for logarithms of norms for the vectors (after orthogonalization but before the nor- 
malization) [21]. From the computations (10 samples, each of 5 ■ 10^ iterations of the 
Poincare map) we obtained the Lyapunov exponents 

Li = 0.6832 ± 0.0007, L2 = -2.6022 ± 0.0036, , . 

L3 = -4.6054 ± 0.0028, L4 = -6.5381 ± 0.0078. ^ ^ 

^For a search of the maximum in a space of three variables at fixed R, we used the Newton method. 

^Restricting the domain of verification of the condition for expanding cones by the set T^(J7) one 
can improve essentially the estimate for maximum allowable 7. In accordance with our computations, 
the inclusion conditions for expanding and contracting cones inside the domain T^([/) are valid yet at 
7^ « 1.5. 
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Presence of the positive exponent Li indicates chaos. (It is close to In 2 = 0.6931 because 
of apphcabihty of the approximate BernouUi map 0.) 

Figure 6 shows portraits of the attractor on a plane of variables of the first oscillator. 
The panel (a) depicts projection of the attractor from the 5D extended phase space on 
the plane of original variables [x, u). The attractor is shown in gray scales (the darkness 
reflects a relative duration of residence inside a given pixel). Black dots relate to the 
Poincare cross-section, the instants tn = nT. The panel (b) shows the attractor in the 
Poincare cross-section on a plane of the redefined coordinates {xq, xi) (see (jl))). Note an 
evident visual similarity with the Smale - Williams attractor as depicted in textbooks. 
The transverse Cantor-like structure is illustrated separately on the panels (c) and (d) 
by magnified fragments of the previous picture. For quantitative characterization of the 
fractal structure in the Poincare cross-section, we have estimated the correlation dimension 
by means of the algorithm of Grassberger and Procaccia. Using a 4-component time 
series x„ = x(t„) obtained from numerical iterations of the Poincare map for n=l-rM, 
M=40000, we get D=1.2516±0.0018 (as a result of averaging over 10 samples). The 
dimension estimated from Luapunov exponents with the Kaplan - Yorke formula is D ~ 
1.263. 

From the point of view of theoretical analysis of the hyperbolic attractors, one of 
the principal features is that intersections of local stable and unstable manifolds if occur 
must be transversal. In computations, to determine the local manifolds with appropriate 
accuracy one can use the following scheme. Let us have three points on the attractor 
obtained one from another by A^-fold application of the Poincare map: ^ = 
T^(xyi) —>■ X(7 = T^(xs), where is a sufficiently large integer. To obtain the ID 
unstable manifold at B, we consider an ensemble of initial conditions close to A and 
parametrized by Aip, a small deflection of the angle variable, of order L^^: xq = r"^ sin?/',, 
Xi = r^cos?/', X2 = x^t x^ = xf, = a/ {xq^ + {x^y, ip = ip"^ + Aip, ip"^ = arg(x^ + 
iXq). After N iterations of the map T, the points take up positions along the unstable 
manifold F^. To obtain the 3D stable manifold at B we set initial conditions for the 
Poincare map close to B: Xq = {r^ + Ar)sin-?/'0) = {f^^ + Ar) cosipQ, X2 = + 
Ax2, = xf + Axs, where = a/ {x^y + (xf )^. Fixing three values (Ar, Ax2, Axs), 
which parametrize the manifold, we take as initial guess ipo = ip^ = arg(xf + iXq) and 
perform N iterations of the map. Then, we get a discrepancy ip^ — ip'" , ipo = tp^ = 
arg(xf -|- IXq), correct the initial angle variable, ipQ = ipo + {ip'^ — '?/'Ar)/2^, and repeat 
the procedure, until the error will be less than a given small value. 

A graphic representation of the manifolds is not trivial because the phase space is four- 
dimensional. Let us use a plane of variables {xq, xi) relating to the first oscillator. The 
ID unstable manifold we show simply as a projection onto this plane. For representation 
of a three-dimensional stable manifold we will use a curve of intersection of the manifold 
with a two-dimensional plane {x2 = , X3 = } projected onto the plane (xq, Xi). 
Practically, a sufficient accuracy for coordinates of points on the manifolds is reached, 
say, at A^ ~ 10. The disposition of the local manifolds revealed from the computations is 
illustrated in Fig. 7. The invariant set that consists of the unstable manifolds coincides 
with the attractor itself. It is enclosed in the toroidal absorbing domain going turn by 
turn around "the hole of the doughnut". On the other hand, the local stable manifolds 
are posed across the "tube" that forms a surface of the toroid. In the two-dimensional 
diagram the stable manifolds look like "speaks of a wheel" . Due to such mutual location. 
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the stable and unstable manifolds can intersect only transversally, and no tangencies do 
occur. 

As stated in this article, in the four-dimensional phase space of Poincare map for 
the system of two non-autonomous coupled van der Pol oscillators there exists a toroidal 
absorbing domain, containing a uniformly hyperbolic attractor. This conclusion is based 
on computer verification of conditions formulated in terms of appropriate inclusion of 
expanding and contracting cones defined in the tangent vector spaces associated with 
the points of the absorbing domain. Hence, our model delivers a long-time expected 
example of a simple physically realistic system with hyperbolic attractor. With this 
example, it will be possible to construct other models with hyperbolic chaos, exploiting 
structural stability intrinsic to the hyperbolic attractors. In fact, a physical experiment 
demonstrating attractor of this type has been performed already on a basis of coupled 
electronic oscillators [16]. In applications, the systems with hyperbolic chaos may be of 
special interest because of their robustness (structural stabihty). An interesting, and now 
a substantial direction is constructing chains, lattices, networks on a base of elements 
with hyperbolic chaos [22]. Models of this class may be of interest for understanding deep 
and fundamental questions, like the problem of turbulence. 

This research was supported by RFBR grant No 06-02-16619. 
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Figure 1: Graphical evidence that the domain U defined by (jH} is absorbing. For initial 
points distributed over a border of U, the resulting data from numerical solution of the 
differential equations over a period T plotted in coordinates ((Tj) fit inside the unit circle 
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Figure 2: The toroidal absorbing domain U and its image T(f/) shown in a 3D projection. 
Variables xq, xi are plotted along the axes in the horizontal plane, and X2 along the vertical 
axis. The fourth variable x^ corresponds to direction of the projecting. 




(a) (b) 

Figure 3: A plot of the function rmax + P = f{R,<^,'4^,(^) at i? = 1 and cj) = 1.25665 
(a) and a diagram illustrating mutual disposition of the 3D cross-sections of the cones 
DTx(5'^) and S'^,. at the point of global maximum of rmax + P- 
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Figure 4: Value of the global maximum of the function / = rmax + p on a hypersurface 
fl26|) versus parameter R (a) and the global maximum value over the whole absorbing 
domain U in dependence on parameter 7. 
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Figure 5: Typical patterns of time dependences for the variables x and y obtained from 
numerical solution of Eqs. Q at T = 6, A = 5, e = 0.5. Panel (a) presents a single 
sample, and panel (b) shows five superimposed samples of the same signal on successive 
time intervals. 
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Figure 6: Portraits of the attractor on a plane of variables of the first oscillator: (a) 
- projection of the attractor from the 5D extended phase space on the plane of original 
variables (x, n); (b) - the attractor in the Poincare cross-section on a plane of the redefined 
coordinates (jH); (c), (d) - details of the Cantor-like transverse structure. 




Figure 7: A diagram on the plane of variables (xq, x\) illustrating mutual location of local 
unstable (u) and stable (s) manifolds for a set of points on the attractor in the Poincare 
cross-section. The gray ring-shaped area depicts a projection of the absorbing domain \J . 
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